arXiv:1506.01539v2 [nlin.AO] 24Jul2015 


APS/123-QED 


Heterogeneity induces emergent functional networks for synchronization 

Francesco Scafuti/ Takaaki Aoki,^ and Mario di Bernardo^’^’0 
^Department of Electrical Engineering and Information Technology, University of Naples Federico II, Naples, Italy 
^Faculty of Education, Kagawa University, Takamatsu, Japan 
^Department of Engineering Mathematics and Bristol Centre for Complexity Sciences, University of Bristol, Bristol, UK 

(Dated: July 27, 2015) 

We study the evolution of heterogeneous networks of oscillators subject to a state-dependent 
interconnection rule. We find that heterogeneity in the node dynamics is key in organizing the 
architecture of the functional emerging networks. We demonstrate that increasing heterogeneity 
among the nodes in state-dependent networks of phase oscillators causes a differentiation in the 
activation probabilities of the links. This, in turn, yields the formation of hubs associated to nodes 
with larger distances from the average frequency of the ensemble. Our generic local evolutionary 
strategy can be used to solve a wide range of synchronization and control problems. 
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I. INTRODUCTION 

Evolution is a fundamental force driving the organi¬ 
zation and structure of natural systems. It is based on 
two key ingredients: variation and natural selection [ij. 
The first ensures the necessary mutation and recombina¬ 
tion generating new species while the second determines 
the survival of the fittest to perform a certain function. 
Networks in Nature have been subject to the same pow¬ 
erful mechanisms that ultimately determined their struc¬ 
ture, properties and functionality. The resulting net¬ 
works have heterogenous topological structures, which 
researchers have been interested in together with their 
effects on dynamical processes Q . Examples include epi¬ 
demic spreading, opinion formation, and synchronization 
ii- Often there is also heterogeneity in the nodes of 
a network. For example, in social networks, individuals 
have different personalities, which will have great impacts 
on their social relationships; or, in manufacturing, indus¬ 
trial products are slightly different from one other, affect¬ 
ing their impact and market shares. The relationship be¬ 
tween the heterogeneity of the nodes and the structural 
properties of a network is little understood, particularly 
when the network evolution is state-dependent. 

Here we suggest that heterogeneity in the nodes is a 
driving force behind the evolution of the network struc¬ 
ture that determines its properties and function. To test 
this ansatz we take as a representative example the prob¬ 
lem of evolving the network structure to achieve synchro¬ 
nization of coupled oscillators. This is one of the best 
understood and most widely studied type of collective 
behavior on networks 

So far, optimal network structures for synchronization 
have been studied mainly by using Monte Carlo meth¬ 
ods [Ml or gradient-based learning strategies [3 [ 3 . 
These are based on the use of some objective function 
for synchronization (as for example the order parameter) 
which is used to find the optimal network whose struc¬ 
tural properties are then surveyed. The Monte Carlo 
approach is a generic and powerful strategy but it is typ¬ 
ically time-consuming, and increasingly cumbersome to 


apply to large-scale networks. Gradient-based methods 
assume some constraints to derive the evolution rule of 
the coupling strengths and the rules are often not local, 
in the sense that some global information on the entire 
network is used. Also, it has been shown that adaptive 
networks can yield the emergence of modular and scale- 
free structures, while enhancing synchronization [3|- 

In this paper, we propose the use of an evolutionary 
strategy to find a functional structure for synchronization 
in a network of heterogeneous oscillators. In so doing we 
will show that heterogeneity in the nodes is instrumental 
in determining the properties of the resulting network. 
The goal of the strategy is to identify, over all possible 
unweighted network configurations, the structure with 
a minimal number of links, which guarantees frequency 
synchronization of its nodes. While the fundamental aim 
of our study is similar to that of the literature [Ml, 
the approach we propose is completely different. Indeed, 
our strategy uses adaptive schemes which are completely 
local and do not rely on any global synchronization mea¬ 
sure. Moreover such schemes are deployed in a novel 
evolutionary manner. 

II. PROBLEM STATEMENT 

We start by considering a network of general nonlinear 
coupled oscillators 
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where Xn € M.P is the p-dimensional state of the n-th 
oscillator, fn denotes its dynamics (note that oscillators 
can be slightly different from each other due to both pa¬ 
rameters and model mismatches), g is a generic coupling 
function and knm are time-varying coupling gains deter¬ 
mining the strength of the coupling between neighboring 
oscillators. 

We model the evolutionary pressures to reach synchro¬ 
nization by considering state-dependent second-order 
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FIG. 1: Schematic description of the evolutionary edge¬ 
snapping strategy. Step 1 (variation): computation of link 
activation probabilities by running the edge-snapping strat¬ 
egy from many different random initial conditions. Step 2 
(selection): selection of those links whose activation probabil¬ 
ity is above some threshold value p*. 


nonlinear dynamics for the gains dependent upon a dou¬ 
ble well potential V{x) = bx'^(x—l)^. The gain dynamics 
are given by 

knm + d knm + - Xn\\), (2) 

in which h(\\xm — a;„||) is a generic increasing function 
such that h{0) = 0. Note that this is a very general adap¬ 
tive network equation relying on a decentralized, local, 
state-dependent interconnection rule. This system can 
be systematically reduced, under a standard technique 
0 , to the network of adaptively coupled phase oscilla¬ 
tors: 
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in which On is the phase of the n-th generic oscillator, 
r(0Tn — On) is a generic 27r-periodic function. We set the 
overall coupling strength K to a unitary value, since it 
can be absorbed into a parameter defining the hetero¬ 
geneity of the natural frequencies by rescaling time, i.e. 
by setting r = Kt. In this paper we analyze, for the sake 
of clarity, the simplest case 
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The effectiveness of edge snapping strategies to achieve 
synchronization has been discussed in |16| and further 
details are given in Appendix A. Under such a forcing 
the dynamics of knm (starting from zero initial condi¬ 
tions knm{0) = 0 and knm{0) = 0); will either converge 
towards 0 (link is not present) or towards 1 (link is acti¬ 
vated) 

The differences in the natural frequencies of the oscilla¬ 
tors originate from the heterogeneity of the node dynam¬ 
ics fn in weakly coupled nonlinear oscillators @ . In what 
follows, these natural frequencies are selected determinis¬ 
tically from a Gaussian distribution with zero mean and 
standard deviation equal to a. Therefore, the parameter 


a can be used to “tune” the level of heterogeneity among 
nodes. 

We note here that when the number of nodes is not 
so large, such as iV = 6 or 7, the natural frequencies 
sampled from a distribution can be biased. To avoid the 
effect of the biased sampling, we deterministically select 
the natural frequencies of the oscillators, similarly to [m, 
as the A^-tuple satisfying the constraints: 

g{uj)duj = (* = 1) 

g{uj)duj = {i = 2,...,N) 



where g{co) is the probability density function of a given 
distribution. It should be noted that for a large network, 
we performed our simulation taking the natural frequen¬ 
cies randomly from a distribution and the obtained re¬ 
sults are qualitatively the same. 

Next, we investigate how the evolution of the network 
is affected by tuning the heterogeneity in the nodes. To 
this aim we use the edge snapping strategy described 
above in a novel evolutionary manner (see Fig. [T]) as 
explained in the next section. 


III. EVOLUTIONARY EDGE-SNAPPING 

The evolutionary Edge-Snapping technique is based on 
two fundamental steps: one implementing the variation 
ingredient of evolution, the other its selection mechanism. 

To implement the variation ingredient of evolution, a 
set of unweighted networks is generated using equations 
([3]) and ([4]) starting the process from different sets of ini¬ 
tial conditions. We consider a set of ns initial conditions 
randomly selected using a Latin Hypercube strategy da 
in the range 0n{0) G [0, 27r[, n = 1,2,..., V. To obtain 
the “fitness” of each link, we next compute the prob¬ 
ability pij of each link being activated as the fraction 
between the number of generated networks where that 
link is present, say riij, and the total number of trials, 
e.g. Pij = nijjns- This yields a stochastic N x N matrix 
P whose elements are the probabilities of activation of 
every possible link among nodes. 

The selection rule is obtained by selecting only those 
links whose activation probability is above a certain criti¬ 
cal threshold value p*, i.e. such that pij > p*. We choose 
p* so as to guarantee that the resulting network is con¬ 
nected and has the smallest number of links. We shall 
term such a network as the minimal edge-snapping (ES) 
network. 

The variation step of our evolutionary strategy relies 
on the generation of a set of ns unweighted network using 
equations m and dl]) and starting the process from a dif¬ 
ferent set of initial condition. With the aim of choosing 
a reasonable value for the number of trials ns, we plot 
in Fig. UKa) the standard deviation of the link activation 
probability pij as a function of ns. As can be noted, the 
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FIG. 2: (a) Standard deviation of the link activation probability pij as a function of the number of trials ns- (b) Order 
parameter and relative number of links of the minimal ES structure as a function of the number of trials ns- 


differentiation in the Pij is quite constant as ns varies 
from 100 to 1000. Thus we select ns = 100 in all of 
our simulations. Indeed this guarantees a good degree of 
variation with the least computational cost. Finally, Fig. 
mb) confirms that the dynamical and structural proper¬ 
ties of the emerging ES minimal structure do not show 
significant fluctuations when the value of ns is increased. 


Note that the state space of the initial phases of 
many oscillators is a high-dimensional space (i.e. the 
aggregate N-dimensional state space obtained collect¬ 
ing the phase of each oscillator in the stack vector 
0 = [ 6 * 1 , 6 * 2 ,..., 6 *Ar]). To obtain effective samplings from 
that space, we adopted a Latin Hypercube Sampling 
(LHS) strategy first proposed in [l^. LHS is a statis¬ 
tical method for generating a sample of plausible collec¬ 
tions of parameter values from a multidimensional dis¬ 
tribution. Specifically, let X denote a N variate ran¬ 
dom variable with probability density function f{x) for 
X € S- Then the range space of each of the N com¬ 
ponents of X is partitioned in ns disjoint intervals Si 
of size Pi = P{X G Si) = 1/ns- Taking the Cartesian 
product of these intervals yields n^ cells each of proba¬ 
bility size ng^ - Each cell can be labeled by a set of N 
coordinates ■ ■ • jiTT-iN) where is the 

interval number of component Xj represented in cell i- 
A LHS is obtained from a random selection of the cells 
TOi,..., mns^ with the condition that for each j the set 
is a permutation of integers 1 , 2 ,..., ns- As a 
result, one random observation is made in each cell. The 
main advantage of the LHS strategy is that it does not 
require more samples for more dimension of the range 
space S- This is the main reason why we use LHS in our 
method. 


To measure the synchronization performance of a ES 
network, we consider an ensemble of phase oscillators 
connected by that network and evaluate Kuramoto or¬ 
der parameter as i?e®^ = ■ 





FIG. 3: (a) Link activation probabilities pij in the case of 
N — 6 generated by the variation stage of the evolutionary 
ES strategy; (b) Selection of the threshold probability value p: 
order parameter R, relative number of links M- The arrow on 
the x-axis indicates the critical threshold p* which gives the 
minimal ES network; (c) Minimal Edge-Snapping Network; 
(d) Optimal network maximizing R obtained by Exhaustive 
search and a Monte Carlo based method. 


IV. EMERGENCE OF MINIMAL NETWORKS 

We first test our strategy by applying it to a small size 
network with N = 6 and cr = 0.3 (Fig. [3]). We obtain 
the P matrix visualized in Fig. 2(a). In Fig. 2(b), as the 
threshold value p is increased, the number of edges, M, 
rapidly decreases while the value of the order parameter 
R remains near unity. 

In the hgure, the normalized number of edges, which 
is divided by maximum links between N nodes, i.e. 
M = is plotted. Also above a certain threshold 

the network becomes disconnected. Therefore we choose 
p* = 0.57 obtaining the minimal ES network depicted 
in Fig. 2(c) which is characterized by M = 7 edges 
and R = 0.96. We compare the minimal ES structure 
with the optimal network structure shown in Fig. 2(d) 
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FIG. 4: Heterogeneity induces functional structural properties of the network. P matrix as a function of the heterogeneity 
parameter a when N — 20. 




(f) 



FIG. 5: Structural properties of the emergent minimal ES network with N = 100 and a = 0.2. (a) Standard deviation of the 
link activation probabilities pij as a function of a. (b) Maximum (red dashed line) and minimum (black solid line) value of 
Node Degree ki as a function of a. (c) Activation probability of each link against the value of the difference between the natural 
frequencies of the oscillators at the endpoints, (d) Node degree ki vs. uii. (e) Order parameter R (red solid line) and relative 
number of links M (blue solid line) of the ES network as a function of the threshold probability value p. For comparison, 
the value R is depicted for an all-to-all network (purple dashed line) and for randomly generated networks (blue dot-dashed 
line) with the same number of links. The arrow on x-axis represents the threshold p* to give the minimal ES network, (f) 
Order parameter R of the phase oscillators interconnected by the minimal ES network when the overall coupling strength K is 
increased (red solid) and decreased (blue dashed). We set N = 300 and cr = 0.2. 


obtained from an exhaustive search and a Monte Carlo 
based method [T^ maximizing the value of R with the 
constraint that the total number of edges M is equal to 
7. We notice that the two networks share the same links. 

Next, we study how heterogeneity induces functional 
structural properties of the network. Figured shows the 
P matrix as a function of the heterogeneity parameter a 
when N = 20. We see that as a is increased a differentia¬ 
tion becomes more and more apparent in the distribution 
of the link activation probabilities pij with edges between 


oscillators with relatively different frequencies becoming 
more likely to occur in the minimal ES structure. 

Fig. [Sja) shows the standard deviation of the link acti¬ 
vation probabilities pij as a linear function of cr in a larger 
network of = 100 oscillators. The structural proper¬ 
ties of the emerging network are therefore induced by the 
node heterogeneity. This is confirmed in Fig. [SKb) where 
the maximum and minimum values of the node degree ki, 
corresponding to each minimal ES network, is plotted as 
a function of a. The behaviors of the maximum value of 
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ki (red dashed line) and the minimum of ki (black solid 
line) show an abrupt transition when passing from a = 0 
to cr > 0. This suggests that the differentiation in the 
degree distribution of the minimal ES network becomes 
remarkable when heterogeneity in the nodes is increased 
from zero (identical oscillators) to a value greater than 
zero (non-identical oscillators). 

The structural properties of the emergent minimal ES 
network are highlighted in Fig. [5jc)-(f) for a network 
of highly heterogeneous N = 100 oscillators (tr = 0.2). 
The activation probability of each link is plotted in Fig. 
[SKc) against the value of the difference between the natu¬ 
ral frequencies of the oscillators at the endpoints. Links 
connecting more distant nodes tend to be activated with 
a higher likelihood confirming that differentiation among 
links is induced by heterogeneity in the nodes. Also, as 
shown in Fig. [SKd), hubs tend to be associated with oscil¬ 
lators whose frequencies are farther away from the aver¬ 
age. The functional advantage of the emerging network 
is shown in Fig. [5{e). Indeed, we observe that the or¬ 
der parameter of the minimal ES structure is close to its 
maximal value for an all-to-all network of the same size, 
even if the number of links in the minimal ES network is 
remarkably lower than that in an all-to-all configuration. 
For the sake of comparison, the values of R for a ran¬ 
domly generated network of the same number of edges is 
also depicted in Fig. EKe). The sudden dip of R is due 
to the graph becoming disconnected beyond that critical 
value of the threshold p*. 

Notice that, as shown in Fig. EKf), as the coupling 
strength K is varied, the order parameter R of the phase 
oscillators interconnected by the minimal ES network ex¬ 
hibits a sudden hysteretic change, associated to a dis¬ 
continuous phase transition, whereas the system with a 
unimodal frequenw distribution undergoes a continuous 
phase transition [^. This discontinuous phase transi¬ 
tion, also known as “explosi ve sy nchronisation”, has been 
studied in the literature [isl - l^ , also in the case of adap¬ 
tive networks [12, [12| , revealing that the correlation be¬ 
tween natural frequencies and the node degree, as shown 
in Fig. 4(d), can induce this phenomenon. Here, we 
wish to emphasise that the proposed evolutionary strat¬ 
egy, which functionally organizes the network structure 
for synchronization, changes the type of phase transition 
that would be generically observed otherwise, inducing 
explosive synchronisation. 

Our results clearly show the role of node heterogene¬ 
ity in inducing functional structures using an evolution¬ 
ary strategy for network synchronization. In particular, 
differences in the node dynamics do influence the evolu¬ 
tion of the network determining a differentiation in the 
link activation probabilities that is instrumental to ob¬ 
tain minimal structures with relatively high values of the 
order parameter. Also, hubs tend to emerge there where 
the distance from the average natural frequency is high¬ 
est. Further simulations also confirmed that a similar 
structure of the emergent network can be induced by us¬ 
ing a power-law rather than a normal distribution when 


selecting the heterogeneous natural frequencies of the os¬ 
cillators (data not shown). 

It is notable that the presence of hubs seems to 
characterize the emergent networks for synchronization 
when the nodes are heterogeneous as opposed to more 
homogenous structures, such as entangled networks, 
which have been suggested to be optimal structures in 
the homogeneous case Q. This is also confirmed in 
the case of Monte Carlo based optimal networks in m 
where the presence of links between nodes with more 
distant frequencies is shown to be more likely and in 
the recent paper [l^ based on the use of gradient-based 
methods. Here we obtain a further confirmation of 
these observations but via a generic local evolutionary 
strategy that is state-dependent and can be applied to 
a wider range of network synchronization and control 
problems. 


V. CONCLUSIONS 

Our results suggest that heterogeneity is the driving 
force determining the evolution of state-dependent func¬ 
tional networks. This can explain the structural proper¬ 
ties detected in natural networks such as neural intercon¬ 
nections in the brain, gene regulatory networks or eco¬ 
logical networks where the states of the nodes t ypically 
affects the evolution of their interconnections [23 - 1^ . 
It can also be used in Dynamical Systems and Control 
theory to design state-dependent evolutionary strategies 
able to induce a desired collective behavior in a network 
of interest. 
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Appendix A: DESCRIPTION OF THE EDGE 
SNAPPING METHOD 


Edge Snapping [T^ is an adaptive strategy for the evo¬ 
lution of an unweighted network. Time-varying coupling 
gains knm are assigned to all pair of nodes n and m, 
with a second order dynamics affected by a double well 
potential V(knm) = bk^^{knm — 1)^ defined as: 


+ d knm + 


dVjknm) 

dknm 


= tl(\\Xm - Xr, 


where d is a damping coefficient, Xn and are the 
states of the nodes at the endpoints of the edge {n,m). 
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FIG. 6: Edges’ evolution according to the edge snapping mechanism. 



FIG. 7: Schematic illustration of how NetEvo works. 


The driving force h{\\xm — a;n||) is a generic increasing 
function such that h(0) = 0. 

The gains’ dynamics mimics the damped motion of a 
particle in a one-dimensional space subject to a double¬ 
well potential as schematically outlined in Fig. [51 In¬ 
deed, in Fig. inija) the initial forcing is strong enough to 
drive the mass particle from the equilibrium at 0 (edge 
turned off) to the well associated to the equilibrium at 
1 (edge activated). On the contrary, in Fig. [Bib) the 


forcing input to the edge snapping dynamics is not able 
to move the particle from the equilibrium at the origin. 
As a result of these dynamics, each coupling gains knm 
converges to either one of the equilibrium points, 0 or 1. 

Note that the dynamics of the gains knm is interdepen¬ 
dent on the dynamics of the nodes’ state (in this let¬ 
ter, the dynamics of the states is given by a coupled oscil¬ 
lator dynamics among nodes). The resulting unweighted 
network is the outcome of the co-evolving dynamics of 
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the nodes and the state-dependent network. This strat¬ 
egy is based on a distributed adaptive nonlinear approach 
and is therefore a generic decentralized approach relying 
only on a nonlinear potential to drive edge adaptation as 
explained in [l^ . 

In addition, this edge-snapping strategy has two pa¬ 
rameters b and d, which can be used to control network 
evolution. Indeed the barrier of the potential between 
the two wells, acts as a constraint. As explained above, 
if the driving force is not strong enough, the edge, after 
a transient, will remain in the well corresponding to the 
absence of link. The height of the barrier can be tuned 
varying the parameter b in the expression of the potential 
V. The higher the barrier b, the stronger the constraint. 

Appendix B: NETEVO 

NetEvo is a computational framework designed to help 
understand the evolution of dynamical complex networks 
[ 2 ^. It provides flexible tools for the simulation of dy¬ 
namical processes on networks and methods for the evo¬ 
lution of underlying topological structures. To bring to¬ 
gether simulation and evolution in a coherent way, the 
framework uses the idea of a supervisor, illustrated in 
Fig. [71 Evolution of the system is performed by the su¬ 
pervisor which can be viewed as a form of optimiser. This 


takes as input an initial topology, simulated output from 
the system and user defined constraints, and aims to re¬ 
turn an optimal or enhanced topology. Changes to the 
system are assessed by using the performance measure -R 
(the opposite of the order parameter), with smaller val¬ 
ues representing an improved performance. By default, 
NetEvo provides a supervisor that uses a Simulated An¬ 
nealing meta-heuristic to search for near optimal configu¬ 
rations. This method has been shown to perform well for 
a wide range of problems with an unknown prior struc¬ 
ture. 

We tuned NetEvo to hnd an optimal structure, given 
an initial condition (the same used in the procedure 
for finding the minimal structure). Simulated annealing 
tends to avoid local minima (or maxima), so we could 
start the optimization from any random connected net¬ 
work structure. However, we decided to start ’’near” the 
minimal structure, to facilitate the optimization (by near, 
we mean a structure obtained from the minimal struc¬ 
ture, after rewiring about 10% of its edges). 

We note here, that it was necessary to run NetEvo sev¬ 
eral times (i.e. n^E = 10 times), because of local max¬ 
ima traps that the algorithm could not avoid. Finally, 
the optimal (or sub-optimal) structure is selected as the 
network that maximize R (starting from Oq), among each 
of the npiE results. 
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